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ABSTRACT 

A  computer  simulation  model  of  the  phosphate,  phytoplank- 
ton  and  zooplankton  dynamics  in  Monterey  Bay  was  examined  and 
modified.   The  model  is  driven  by  four  forcing  functions  ex- 
pressed as  annual  cycles  of  upwelling  velocity,  incident 
solar  radiation,  mixed  layer  depth,  and  mixed  layer  tempera- 
ture.  An  alternate  upwelling  index  was  developed  based  on 
the  local  wind  field.   A  revised  radiation  index  is  employed 
based  on  the  generation  of  both  advective  fog  and  low  stratus 
cloud  cover  common  during  upwelling  on  the  California  coast. 
Analysis  of  the  model's  response  to  sinking  and  advection  of 
phytoplankton  was  examined.   The  importance  of  seasonal  in- 
creases in  predators  was  introduced  as  a  controlling  factor 
in  the  seasonal  growth  of  zooplankton.   The  model  is  able  to 
predict  the  seasonal  trends  of  phosphate,  phytoplankton,  and 
zooplankton  throughout  the  year. 
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I.   INTRODUCTION 

A.   PURPOSE 

In  the  interests  of  assessing  the  biomass  of  oceanic 
areas,  various  mathematical  models  have  been  designed  to 
predict  the  dynamic  response  of  ecosystems  to  change  in  the 
physical  and  chemical  environment.   Such  models  necessarily 
rely  on  accurate  characterization  of  the  forcing  functions 
and  accurate  representation  of  the  system  with  equations 
which  are  consistent  with  the  specific  region  under  study. 

Recent  attention  paid  to  modeling  systems  in  upwelling 
regions,  Coastal  Upwelling  Ecosystems  Analysis  (CUEA),has 
brought  about  developments  in  both  these  areas.   Adequate 
data  of  a  physical  and  chemical  nature  is  becoming  available 
and  significant  progress  in  the  refinement  of  governing  equa- 
tions is  evident  (Walsh,  1973). 

The  research  presented  here  was  aimed  at  creating  a  simu- 
lation model  to  describe  the  seasonal  plankton  dynamics  in 
Monterey  Bay,  California,  as  known  from  the  best  available 
data.   An  existing  simulation  model  (Pearson,  1975)  was 
evaluated  in  an  effort  to  make  refinements  and  to  judge  the 
applicability  of  the  time  simulation  technique  to  the 
Monterey  region,  characterized  by  seasonal  upwelling. 

3.   MODELING  PHILOSOPHY  IN  THE  ECOLOGICAL  CONTEXT 

The  merit  of  mathematically  modeling  systems  of  a  purely 
physical  nature  in  which  parameters,  initial  conditions, 
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boundary  conditions  and  the  variable  relationships  are  readi- 
ly specified  is  unique.   The  value  of  current  models  of  a 
biological  nature  requires  some  comment. 

According  to  Odum  (1971)  there  are  a  number  of  reasons 
behind  constructing  any  mathematical  model.   Prediction  of 
future  states  of  the  variables  involved  is  an  end  in  itself. 
Use  of  the  prediction  accuracy  to  focus  data  acquisition 
efforts  or  direct  attention  to  deficiencies  in  concepts  is 
likewise  valid.   If  the  model  is  real  in  the  sense  that  it 
attempts  to  predict  changes  in  the  state  variables  based  on 
a  framework  drawn  from  research  in  the  real  world,  then 
failures  may  be  useful  in  delineating  weaknesses  in  the 
governing  equations. 

Construction  of  a  "real"  model  is  at  best  difficult.   A 
perfect  model  of  an  ecosystem  assumes  that  "all  states  and 
rates  of  change  of  the  variables  are  known  at  all  times" 
(Walsh,  1973).   Given  the  vastness  of  complex  interactions 
in  even  the  simplest  of  biological  systems,  it  may  be  stated 
that  "no  perfect  representation  of  the  real  world  exists 
except  the  real  world  itself"  (Walsh,  1973). 

It  follows  then  that  certain  assumptions,  logical  state- 
ments and  approximations  of  real  dynamics  are  necessary  to 
allow  the  model  to  operate.   The  resulting  model  may  bear 
little  correlation  to  the  actual  physical,  chemical  and 
biological  situation  at  hand.   Again,  the  argument  for  con- 
structing the  model  must  be  considered.   For  purposes  of 
prediction,  a  model  which  makes  a  logical  statement  relating 
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two  state  variables  in  the  attempt  to  specify  a  relationship 
which  may  be  unknown  in  part  or  in  whole  has  merit.   The 
fact  that  the  logic  of  the  model  is  a  crude  simplification 
is  incidental  if  the  prediction  capability  is  proven. 

All  models  can  be  characterized  on  the  basis  of  realism, 
precision,  and  generality  (Odum,  1971).   It  has  been  suggested 
that  these  three  qualities  are  mutually  exclusive  in  ecologi- 
cal models  (Patten,  1971).   A  fourth  category,  simplicity, 
might  be  considered.   Exacting  detail,  while  lending  to 
realism  and  precision  often  forces  the  model  to  be  specific. 
Only  when  the  model  becomes  an  isomorph,  i.e.,  a  one-to-one 
correspondence  to  the  real  world,  will  generality  be  restored 
(Walsh,  1973). 

The  problems  in  creating  a  realistic  model  have  been  dis- 
cussed briefly.  The  value  of  this  type  of  model  lies  in  re- 
search guidance.   Where  a  specific  question  is  asked  of  an 
ecological  model,  precision  in  prediction  capability  may  be 
enhanced  by  sacrificing  the  other  qualities.   Such  an 
approach  might  be  taken  in  a  fisheries  model  where  a  precise 
output  for  a  limited  region  is  desired  (Odum,  1971).   Simplic- 
ity and  generality  at  the  expense  of  both  precision  and 
realism  may  enhance  understanding  of  broad  ecological  con- 
cepts . 

The  study  presented  herein  was  conducted  with  the  idea 
of  creating  a  model  having  good  predictive  capabilities  for 
conditions  occuring  over  a  long  span  of  time  in  a  limited 
oceanic  region. 
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C.   DESCRIPTION  OF  THE  MODEL 

The  original  ecosystem  model  was  developed  by  Pearson 
(1975)  as  part  of  continuing  research  into  the  correlation 
of  zooplankton  biomass  with  the  chemical  and  acoustic  pro- 
perties of  the  ocean  (Traganza,  1976).   It  is  the  intent  of 
this  research  to  identify  areas  of  weakness  in  the  model 
and  make  the  appropriate  changes  to  the  equations  and  input 
parameters.   These  changes  include  a  refinement  of  the  up- 
welling  index  and  radiation  index  as  well  as  inclusion  of 
advection,  sinking  and  predation  terms  in  the  governing 
equations . 

The  modeling  technique  used  is  a  time  dependent  simula- 
tion solved  digitally  on  an  IBM  360  computer  using  the  IBM 
Continuous  Simulation  Modeling  Program  (CSMP-360).   The 
dynamic  equations  are  written  in  non-linear  differential 
form  and  driven  by  four  exogenous  forcing  functions  expressed 
as  annual  cycles  of  upwelling  velocity,  mixed  layer  depth, 
mixed  layer  temperature,  and  incident  solar  radiation.   The 
output  consists  of  annual  cycles  of  phosphate  concentration 

Cmicrogram-atoms  of  phosphorous  per  liter,  ug-atP/£),  phyto- 

2 
plankton  biomass  (grams  of  carbon  per  square  meter,  gC/m  ) 

2 

and  standing  stock  of  herbivorous  zooplankton  (g  C/m  )  in 

the  mixed  layer. 

The  Pearson  version  (Figure  1)  of  the  model  is  defined 
by  the  following  basic  equations: 

1)  IP-  =  NUT  +  REGEN  -  UPTAK  (1) 

dt 
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LIGHT 


Figure  1.   Pictoral  representation  of  Monterey  Upwelling 
Ecosystem. 
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whe^e  *     dXl 

-^r—   =  time  rate  of  change  of  phosphate  concentra- 
tion (ug-at  P/l)  in  the  mixed  layer 

NUT  =  input  of  phosphate  from  upwelling  and  mixing 
from  below  the  mixed  layer  (yg-at  P/l) 

REGEN  =  phosphate  recycled  as  biological  excretory 
products  within  the  mixed  layer  (yg-at  P/l) 

UPTAK  =  phosphate  depleted  by  phytoplankton  utiliza- 
tion (yg-an  P/l) 

2)  =¥-   =  PROD  -  RESP2  -  GRAZ  (2) 
dt 

dX2 

where:     -7- —  =  time  rate  of  change  of  phytoplankton  biomass 

dt     (gC/m2) 

2 
PROD  =  photosynthetic  production  (gC/m  ) 

2 
RESP2  =  respiration  losses  (gC/m  ) 

2 
GRAZ  =' grazing  losses  (gC/m  ) 

3)  5S-1  =  GRAZ  -  RESP3  -  VOID  -  LOSS  (3) 
dt 

GRAZ  =  input9due  to  ingestion  of  phytoplankton 
CgC/nT) 

2 

RESP3  =  respiration  losses  (gC/m  ) 

2 
VOID  =  excretory  losses  (gC/m  ) 

2 
LOSS  =  stock  depletion  by  mortality  (gC/m  ) 

(Pearson,  1975) 
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II.   BACKGROUND  THEORY 

A.   UPWELLING  THEORY 
1.   Ekman  Dynamics 

Seasonal  upwelling  of  nutrient-rich  deep  water  into 
the  productive  mixed  layer  is  an  important  aspect  of  the 
•Monterey  Bay  ecosystem.   The  principal  source  of  phosphate 
in  the  mixed  layer  is  vertical  advection  which  is  produced 
in  response  to  wind-induced  Ekman  circulation.   Theoretical 
calculations  of  the  offshore  component  of  horizontal  current 
are  based  on  Ekman  pure  drift  theory  in  an  infinitely  deep, 
homogeneous  ocean  (Neumann  and  Pierson,  1966).   The  compo- 
nents of  the  horizontal  current  take  the  form: 

U   =    V.     e-(7T/D)z    C0S[45-U/D)Z]       and  (4a) 

o 

V    =    V      e"(TT/D)Z    SIN[45-(tt/D)Z]  (4b) 

o 

D  =  (A/pu)  SIN  <|>)'s   and  (5a) 

V   =  T  /Aa  (5b) 

o     y 


where 


U  =  velocity  component  in  the  x  direction 

V  =  velocity  component  in  the  y  direction 
(parallel  to  wind) 

V   =  surface  velocity  component  (45   to  right  of 
wind j  Northern  Hemisphere) 

Z  =  depth 

D  =  Ekman  depth  of  frictional  resistance 

A  =  coefficient  of  eddy  viscosity 
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t   =  surface  wind  stress  in  y  direction  (dynes) 
a  =  fp/A 

<p    =  latitude 

3 
p  =  density  of  water  (g/cm  ) 

f  =  Coriolis  parameter  =  2w  SIN  <j> 


The  coefficient  of  eddy  viscosity,  A,  was  assumed 
constant  since  detailed  information  on  the  small  scale  mo- 
tions of  the  surface  layer  was  not  available.   This  assump- 
tion is  generally  made  in  mass  transport  studies  (Neumann, 
1968). 

The  mass  transport  due  to  wind  driven  current  is 
found  by  integrating  the  equations  of  velocity  (U  and  V) 
over  the  depth  of  the  water  column,  as  shown: 


S   =  p    UdZ  =  t  /f  (6a) 

x  y 

and     S   =  p    VdZ  =  0  (6b) 

y 

where : 

S   =  mass  transport  in  the  x  direction 
S   =  mass  transport  in  the  y  direction 

From  the  equations  for  mass  transport,  it  is  observed  that 
the  net  transport  is  90  degrees  to  the  right  of  the  wind  in 
the  northern  hemisphere  and  that  it  is  directly  proportional 
to  the  wind  stress  acting  on  the  sea  surface,  parallel  to 
the  coast. 
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2 .   Wind  Stress  on  the  Sea  Surface 

The  coupling  of  wind  energy  to  the  water  at  the  air- 
sea  interface  is  defined  by  the  conventionally-accepted  form 
of  the  wind  stress  equation : 

T   =  PfCzW^  (Gordon,  1972)        (7) 

2 
where:   T   =  wind  stress  parallel  to  the  coast  (dynes/cm  ) 

p1  =  density  of  air 

C   =  drag  coefficient  at  height  z  (non-dimensional) 

z 

W   =  wind  speed  at  height  z  (cm/sec) 

As  in  other  momentum  exchange  studies ,  the  dynamic 
coefficients  appear  to  present  the  singular  most  difficult 
problem.   Wilson  (1960)  states  that  the  drag  coefficient  for 
air  flow  over  a  water  surface  is  dependent  on  the  wind  speed 
(W) ,  the  height  of  measurement (z ) ,  a  surface  roughness  para- 
meter (z  ),  and  atmospheric  stability  terms.   There  may  be 
additional  dependence  on  oceanic  parameters  of  depth,  fetch 
distance  and  wave  height.   An  average  value  of  the  drag  co- 
efficient for  high  wind  speeds  (i.e.  greater  than  10  m/sec) 

-3 
of  2.37  X  10    is  arrived  at  by  examination  of  research  com- 
pleted through  1959  (Wilson,  1960). 

Work  by  Deacon  (1962)  established  an  empirical  rela- 
tionship for  winds  up  to  13  meters/second  which  gives  a 
linear  dependence  of  C   (at  z  =  10m)  to  the  wind  speed  as 
follows : 

C1Q  =  (1.10  +  0.04  W   )  X  10"3  for  W  in  m/sec       (8) 

(Roll,  1965) 
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Deacon's  research  was  carried  out  using  data  from  ship  obser- 
vations and  coastal  regions. 
3.   Spatial  Extent 

The  offshore  movement  of  water  in  a  region  bounded 
by  a  coastline  causes  replacement  water  to  be  upwelled  from 
below  the  layer  affected  by  the  surface  wind  stress.   Accord- 
ingly, it  is  necessary  to  specify  the  spatial  dimensions  of 
the  region  in  order  to  apply  principles  of  continuity  in  de- 
termining the  vertical  current  speed.   The  significant  para- 
meters involved  in  upwelling  are:  (1)  the  offshore  horizontal 
dimension  of  the  upwelling  region,  which  is  most  directly 
related  to  the  width  of  the  wind  field;  (2)  the  offshore 
component  of  the  Ekman  current  (U) ;  and  (3)  the  depth  of  the 
Ekman  layer  (D) . 

Estimates  of  the  maximum  depth  of  the  source  of  up- 
welled  water  in  a  coastal  region  are  approximately  100  to 
200  meters  (Neuman  and  Pierson,  1966)  based  on  the  slope  of 
the  isotherms.   The  mass  transport  offshore  due  to  wind 
drift  current  has  been  previously  detailed  as  occurring  in 
a  surface  layer  of  depth,  D,  corresponding  to  the  Ekman  depth 
of  frictional  resistance;  regardless  of  the  source  depth  of 
the  upwelled  water,  it  transgresses  the  "boundary"  at  depth, 
D,  before  being  carried  offshore.   The  vertical  dimension  of 
the  region  is  then  readily  specified.   It  is  assumed  that 
the  depth  of  frictional  resistance  (D)  and  the  depth  of  the 
mixed  layer  (Z)  are  approximately  coincident.   The  differ- 
ences between  (D)  and  (Z)  are  important  when  phytoplankton 
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are  advected  out  of  the  Ekman  layer  but  exist  throughout 
the  mixed  layer. 

The  horizontal  extent  of  the  upwelling  region  is 
considerably  more  nebulous  than  the  depth  of  the  mixed  layer. 
Classical  estimates  of  this  dimension  limit  the  zone  to 
probably  no  more  than  10  0  kilometers  in  offshore  width.   It 
may  be  expected  that  the  width  of  the  region  is  dependent 
on  several  factors;  among  them,  the  width  of  the  wind  field, 
variations  and  non-linearities  in  the  energy  exchange  pro- 
cesses due  to  surface  roughness  characteristics  or  stability 
changes  in  both  the  atmosphere  and  the  water  column,  as  well 
as  spatial  and  temporal  oscillations  of  the  significant  wind 
vector.   In  addition,  local  coastal  and  bottom  topography 
may  figure  extensively  in  the  problem  as  will  be  discussed 
later. 

Sverdrup  (1938)  observed  a  relatively  well  defined 
offshore  boundary  to  a  coastal  upwelling  zone,  coincident 
witha  downwind  current  which  was  marked  by  an  intense  verti- 
cal gradient  of  velocity  (Sverdrup  et  al.,  19M-2).   He  further 
observed  an  offshore  migration  of  the  current  band  at  a  rate 
somewhat  less  than  the  speed  of  the  offshore  surface  current 
within  the  upwelling  region.   The  latter  fact  implied  the 
possibility  of  cellular  circulation  patterns  in  the  near  sur- 
face upwelling  zone. 

Hidaka  (195M-)  proposed  a  steady  state  theory  of  up- 
welling in  which  he  defined  a  horizontal  frictional  distance, 
D,  ,  analogous  to  the  depth  of  frictional  resistance 
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appearing  in  Ekman's  work.   Applying  reasonable  values  to 
Hidaka's  expression  confines  the  total  cellular  circulation 
region  to  a  width  of  about  15  kilometers.   The  width  of  the 
region  is  defined  by: 

DR    =    tt(2Ah/2o)    SIN    <J>)2  where:    AR    =    107(gcm_1    sec'1)         (9) 

(Smith,  1968) 
Downward  vertical  currents  occurring  in  the  seaward  half  of 
the  cell  limit  the  upwelling  to  a  width  of  7.5  kilometers. 
The  theoretical  width  is  not  consistent  with  observations  of 
upwelling  at  greater  distances  from  the  coast  (Barnes,  1969). 

Yoshida  (1967)  developed  a  quasi-steady  state  model 
applicable  to  eastern  boundary  current  regions  (Smith,  19  6  8^. 
An  expression  was  derived  for  the  horizontal  dimension  of 
the  coastal  upwelling  region  which  is  given  by: 


L   =  (gH^p/p)   (Hy/Dy)^  for  latitudes  greater  than 
x      x  22  degrees  (10) 


where:    L   =  horizontal  width  of  the  coastal  boundary 
region  (m) 

_2 
g   =  gravitational  acceleration  (m-sec   ) 

H   =  thickness  of  the  upper  layer  (m) 

D   =  H  +  thickness  of  the  lower  layer  (m) 

f   =  Coriolis  parameter  (sec"  ) 

_3 
p   =  density  of  water  (g-cm   ) 

Ap   =  density  difference  between  deep  and  surface 
layers  (g-cm-3) 

y   =  internal  friction  coefficient  (sec   ) 

Y   =  2  X  10"7  sec"1  (Smith,  1968) 
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3.   CURRENTLY  AVAILABLE  UPWELLING  INDEX 

Bakun  (197  6)  has  calculated  coastal  upwelling  indices 
for  the  coast  of  North  America  using  Fleet  Numerical  Weather 
Central  monthly  mean  pressure  fields  to  compute  the  geos- 
trophic  wind.   Analysis  is  done  on  a  63  x  63  point  grid  of 
approximately  200  nautical  mile  mesh  length  (Bakun,  1973). 
Figure  2  shows  Bakun' s  results  for  1974  at  36  degrees  North, 
122  degrees  West,  approximatley  5  0  miles  south  of  Monterey 
Bay.   In  reviewing  several  aspects  of  hydrographic  surveys 
taken  in  Monterey  Bay,  it  was  suspected  that  this  region  is 
characterized  by  local  upwelling  patterns  not  appearing  in 
the  index  computed  by  Bakun.   Further  discussion  follows  in 
the  methods  section. 

C.   EFFECTS  OF  UPWELLING  ON  INCIDENT  RADIATION 

Cold  water  upwelled  from  depth  during  the  spring  and 
summer  months  brings  about  both  advection  fog  and  stratus 
cloud  formation.   The  fog  occurs  as  warm  surface  air  is 
cooled  by  the  cold  seawater  to  its  saturation  point  causing 
the  condensation  of  the  contained  water  vapor.   Stratus 
clouds  occur  below  the  base  of  a  quasi-permanent  atmospheric 
temperature  inversion  (Tont,  197  5). 

The  net  effect  of  both  fog  and  stratus  layers  is  to  re- 
duce the  amount  of  solar  radiation  reaching  the  sea  surface 
during  the  upwelling  period.   Tont  (1975)  obtained  mean 
values  of  solar  radiation  at  the  surface  over  the  period 
1950-1973  at  San  Diego.   A  correlation  study  of  the  annual 
upwelling  index  and  the  amount  of  sunlight  at  the  surface 
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showed  a  maximum  of  about  M-0  per  cent  reduction  in  available 
sunlight  during  the  peak  upwelling  months.   The  result  of 
Tont ' s  research  is  shown  in  Figure  3. 

D.   SINKING  OF  ALGAL  CELLS 

Review  of  sedimentation  rates  of  algal  cells  indicate 
justification  for  including  a  sinking  term  in  the  phytoplank- 
ton  equation.   The  Pearson  (1975)  model  assumes  uniform  algal 
concentation  within  the  mixed  layer.   It  is  reasonable  that 
changes  in  the  concentration  of  phytoplankton  caused  by  ver- 
tical movement  out  of  the  layer  will  occur. 

Parsons  and  Takahashi  (1973)  state  that  the  theoretical 
sinking  rate  of  phytoplankton  can  be  determined  from: 

v  =  ifli  (£l=£)  (ID 

where : 

r  =  radius  of  the  cell 

g  =  gravitational  acceleration 

p'  =  density  of  the  organism 

p  =  density  of  the  medium  (sea  water) 

r\  -    viscosity  of  the  medium 

<J>   =  form  resistance  coefficient  (accounts  for  non- 
r 

spherical  shapes) 

Measured  rates  of  sinking  for  live  algal  cells  vary  be- 
tween zero  and  3  0  meters  per  day  according  to  Parsons  and 
Takahashi  (1973).   These  rates  represent  a  broad  range  of 
soecies  and  sizes. 


24 


3N1HSNPIS     3"l9ISS0d      lN30d3d 


o 

CD 

& 

CO 

CD 
-H 

CO 

4h 

£ 

«d 

v^ 

s_^ 

X 

CD 

CD 

C 

T3 

•H 

C     rC 

■H 

bO 

3 

C 

w 

•H 

H 

CD 

H 

H 

CD 

£1 

3 

•H 

ft  W 

3 

en 
0 

Uh 

ft 

o 

4-1 

01 

0 

c 

• 

(d 

■M 

/— \ 

CD 

C 

LO 

6 

CD 

r- 

0 

CD 

>,  f- 

H 

CD 

,C 

ft 

*N 

-M 

-P 

c  -n 

c 

0 

C 

0 

s 

rd 

H 

• 

00 

CD 

h 

3 

b£ 

•H 

Fm 

O 

o 

o 
o 

i-oas-puj) 

X3QN 

o 
o 


9NI"113Mdn 


25 


E.   ADVECTION 

It  is  known  that  significant  current  velocities  occur 
within  the  mixed  layer  from  Ekman  dynamics.   It  might  be  ex- 
pected then  that  the  concentration  of  phosphate,  phytoplank- 
ton  and  zcoplankton  will  be  affected  by  advection.   O'Brien 
and  Wrobleski  (1972)  have  shown  by  scale  analysis  that  advec- 
tion and  biological  productivity  of  the  Peru  upwelling  eco- 
system are  equally  important. 

The  advective  terms  of  the  material  derivative  of  a 
property  (C)  are  written: 

Advective    change    =u|§+v|y+WI§  (12) 


where:   U,V,W  =  velocity  components  in  x,  y,  and  z  directions 

respectively 

3C   3C   3C 

ttt,  iry,   -yY  -    gradients  of  the  concentration  of  property 

(C)  in  the  x,  y,  and  z  directions 


A  review  of  California  Cooperative  Fisheries  Investiga- 
tion (CALCOFI)  data  for  1974  (CALCOFI,  1974)  shows  that  the 
phosphate  concentration  varies  significantly  in  the  vertical 
direction  but  the  horizontal  variation  is  not  as  significant. 

It  is  assumed  for  the  purposes  of  this  model  that  the  ad- 
vection term  for  phosphate  reduces  to  the  vertical  component: 

8X1 
Advective  change  =  W  js—  (13) 

3X1    3X1    n 
since:        ^-  =  ^-   -  0 

Pearson  (1975)  has  included  this  term  in  the  computer  model 

as:  UPWELL  =  W  DX1  ~  X1  (14) 

Jl  J\J_i 
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where : 

W  =  vertical  speed  (m/day) 

EKL  =  depth  of  the  Ekman  layer  (m) ,  EKL  =  D 

DX1  =  concentration  of  phosphate  below  the  mixed 
layer  (yg-at  P/l) 

XI  =  concentration  of  phosphate  in  the  mixed  layer 
(constant)  (yg-at  P/l) 

DX1  -  XI 

gY? =  vertical  gradient  of  phosphate  (yg-at  P/m  1) 

The  advection  of  zooplankton  in  the  revised  model  is 
discounted  because  it  is  assumed  that  they  have  developed 
mechanisms  which  permit  them  to  remain  in  a  particular 
oceanic  region  by  swimming  or  riding  cellular  circulation 
currents . 

Phytoplankton,  however,  are  affected  by  advection. 
CALCOFI  data  (CALCOFI,  1974)  shows  that  the  significant 
gradients  of  phytoplankton  occur  in  the  vertical  and  off- 
shore directions;  therefore,  the  advection  expression  for 
phytoplankton  is  reduced  to: 

Advective  change  =  U  rrrr—   +  W  -^y—  (15) 


The  horizontal  velocity  is  related  to  the  vertical  veloc- 

dV 

dY 


ity  by  continuity  (when  it  is  assumed  that  -jy-  =  0,  parallel 


to  the  coast)  such  that: 

U  =  -  WL/D  (16) 
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where:      U  =  horizontal  velocity  component  offshore  (m/day) 
W  =  upward  vertical  velocity  component  (m/day) 
L  =  width  of  the  region  (m) 
D  =  Ekman  depth  (m) 

The  phytoplankton  advection  term  can  be  written: 


Advective  change  =  W  [-  |  ||^  +  ||i]  (17) 

The  gradients  of  phytoplankton  can  be  approximated  by: 


9X2     LIM  AX2       ,  3X2  _   LIM  AX2  MQ  ... 

9X~  '  AX+0  AX~     ana  3Z~  "  AZ+0  AZ~  Qiaafet; 


If  it  is  assumed  that  the  concentrations  of  phytoplankton 
below  the  mixed  layer  and  outside  the  modeled  area  are  much 
less  than  the  concentration  in  the  modeled  region,  then  the 
gradient  can  be  approximated  by: 


=  horizontal  gradient;  (19) 


AX 

X2-0 


(20) 


■j-j—  -    vertical  gradient; 

where:  X2  =  phytoplankton  concentration  in  the 

mixed  layer 

AX  =  horizontal  distance  at  which  the  concen- 
tration of  phytoplankton  becomes  negli- 
gible with  respect  to  the  modeled  area 

AZ  =  vertical  distance  at  which  phytoplankton 
concentration  becomes  negligibly  small 

The  phytoplankton  advection  term  can  then  be  written  as : 

Advective  change  to  ,      , 

Phytoplankton  in  the  =  W(X2)[-  ~%   +  aY]  (21) 

modeled  area 
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The  sign  convention  is  established  as  positive  (+)  into  the 
model  area  and  negative  (-)  out. 

F.   PREDATION  PRESSURE 

Under  constant  environmental  conditions ,  the  dynamics  of 
zooplankton  growth  are  a  function  of  the  food  supply  and 
predation  pressure.   Natural  processes  of  population  control 
insure  that  a  given  species  will  not  be  eliminated  by  low 
food  supply  and  high  predation  pressure  or  that  unstable 
growth  will  not  result  from  the  inverse  situation.   Attempts 
to  translate  the  stabilizing  factors  into  mathematical  terms 
often  result  in  simplified  general  relationships  that  do  not 
exhibit  the  flexibilities  inherent  in  the  real  system.   Not- 
withstanding the  limitations,  a  predation  term  is  included 
in  the  simulation  model. 

The  long  term  survivability  of  both  predator  and  prey  is 
keyed  directly  to  maintenance  of  a  balance  in  energy  expendi- 
tures and  gains.   To  date,  laboratory  experiments  to  dupli- 
cate this  balance  have  not  been  entirely  successful.   Prey 
stocks  in  small  laboratory  environments  have  been  artificial- 
ly supported  by  providing  refuge  where  predator  access  to 
the  prey  was  denied  and  by  introducing  additional  prey  to 
the  system  to  replenish  the  supply  (Patten,  1971).   Neverthe- 
less, population  control  by  predation  is  prey-density  depen- 
dent as  a  first  approximation.   For  a  system  of  one  prey 
species  and  a  single  predator  species,  predation  pressure 
may  be  represented  by: 

P  =  k1  D   for  0<D<d1  and 

P  =  k2     for  B>d1  (22) 
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where:     P  =  predation  pressure  or  food  intake  per  unit 
predator 

D  =  prey  density 

d  =  satiation  threshold  of  prey  density 

k,  and  k«  =  constants 

Above  prey  density  values  of  d, ,  predator  satiation  will  oc- 
cur and  pressure  will  level  off  at  some  constant  high  value. 
The  function  is  shown  on  an  arbitrary  scale  in  Figure  4 
according  to  Patten  (1971).   This  simple  relationship  is 
not  entirely  satisfactory  for  the  predator  since  low  food 
supplies  would  result  in  starvation.   Under  such  circum- 
stances, and  given  a  second  prey  species  of  perhaps  less 
palatability ,  the  predator  would  switch  food  intake  to  the 
more  abundant  species  of  prey.   There  are  several  obvious 
implications  to  the  switching  phenomenon: 

1.  predation  on  the  most  abundant  species  assures 
stable  growth  of  that  species , 

2.  switching  from  a  declining  species  serves  to  prevent 
elimination  of  that  species  (Patten,  1971),  and 

3.  elimination  of  the  predator  by  overgrazing  one  avail- 
able food  source  is  avoided. 

The  model  under  consideration  in  this  thesis  is  a  single 
species  approximation  in  which  predation  pressure  is  incor- 
porated with  other  zooplankton  mortality  factors,  i.e. 
morphological  death,  in  the  "LOSS"  term.   This  term  is  given 
by  Pearson  (1975)  as: 

LOSS  =  L  (X3)  (23) 
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PREY     DENSITY 


Figure  U- .   Predation  pressure  as  a  function  of  prey 
density  (zooplankton  biomass)  (after 
Patten,  1971). 
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where:     LOSS  =  fraction  of  the  zooplankton  population  re- 
moved from  the  system  in  a  day's  time 
(g  C/m2  day) 

L  =  constant  rate  of  loss  or  percent  loss  per 
day  (day~l) 

o 
X3  =  zooplankton  "standing  crop"  (g  C/m  ) 


Harriston  concluded  in  1960  that  "herbivores  are  preda- 
tor limited"  (Patten,  1971). 

Conclusions  by  Patten  (1971)  show  that  herbivores  are 
both  food  and  predator  limited.   Therefore,  in  the  modified 
model,  the  LOSS  term  is  reserved  for  variable  predator  limi- 
tation of  herbivorous  zooplankton  and  natural  death  of 
zooplankton  herbivores  is  assumed  to  be  negligible. 
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III.   METHODS 

The  existing  model  appears  weak  in  several  areas.   Orig- 
inally, the  forcing  functions  of  upwelling  and  solar  radia- 
tion which  are  critical  to  the  simulation  were  obtained  from 
average  cycles  for  the  California  coast  (upwelling)  and  the 
latitude  band  between  3  0  and  M-0  degrees  north.   Phytoplankton 
sinking  and  advection  terms  were  lacking  and  the  predation 
pressure  term  in  the  zooplankton  equation  needed  improvement. 

A.   UPWELLING  PROGRAM 

A  computer  program  was  written  to  calculate  the  annual 
upwelling  index  for  Monterey  Bay.   Wind  data  was  obtained 
from  the  Monterey  Peninsula  Airport  observations  during  1974. 
Hourly  observations  were  examined  to  determine  the  signifi- 
cant mean  daily  wind  vector  for  each  day  and  corrections 
were  applied  to  speed  and  direction  when  indicated  by  topo- 
graphic obstructions.   Fortunately,  no  corrections  were 
needed  to  northerly  winds  which  are  the  driving  mechanism 
for  upwelling. 

An  x-y  coordinate  system  was  established  such  that  the 

2 

surface  wind  stress  (x   =  p'C  W  )  could  be  directly  calcu- 

y 

lated  from  the  wind  component  parallel  to  the  coast.   The 
Ekman  mass  transport  (S  ) ,  surface  current  velocity  (V  ) 
and  the  average  velocity  within  the  Ekman  layer  (U)  were 
calculated  by  the  following  relations : 

V   =  x  /Aa  (5b) 

o    y 
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where:      V   =  surface  current  vector 
o 


and 


■»/ 


U  =  ±  /  UdZ  (24) 


where:       U  =  average  velocity  in  x  direction  in  the  layer 
and: 

S   =  x  /f  (6a) 

x     y  v   ' 

where      S   =  mass  transport  in  x  direction  (offshore) 

The  vertical  current  representing  upwelling  is  determined 
per  unit  y  (along  coast)  from  the  continuity  equation, 
(7  -  v)  =  0,  such  that 

0(D)  =  W(L)  (25) 

where:      U  .=  average  horizontal  velocity  in  the  Ekman 

layer  (m/day) 

W  =  upwelling  speed  (m/day) 

D  =  Ekman  depth  (m) 

L  =  horizontal  width  of  the  upwelling  region  (m) 
Zero  current  is  assumed  parallel  to  the  coast. 

The  horizontal  width  of  the  upwelling  region,  (L)  was 
determined  from  Yoshida's  (1967)  equation.   Applying  typical 
values  to  this  equation  yields  a  horizontal  width  of  3  X 
10   m. 

Principles  of  continuity  were  used  to  arrive  at  the  mean 
daily  upwelling  current  values  and  these  figures  were  aver- 
aged over  a  3  0-day  period. 
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B.   RADIATION  INDEX 

A  revised  incident  solar  radiation  index  which  takes  in- 
to account  the  reduction  of  sunlight  by  local  fog  and  stratus 
cloud  cover  was  developed.   Correlation  of  the  measurements 
of  surface  irradiance  and  the  strength  of  upwelling  (as  in- 
dicated by  the  upwelling  index)  by  Tont  (1975)  is  repre- 
sented by  Figure  5.   This  figure  is  derived  from  the  results 
of  Tont ' s  study  as  shown  in  Figure  3  by  plotting  the  percent 
of  possible  sunshine  (Y-axis)  as  a  function  of  the  upwelling 
index  (X-axis)  which  has  been  scaled  to  a  maximum  value  of 
one.   The  resulting  curves  (A)  and  (B)  show  the  conditions 
occurring  after, and  prior  to,  the  upwelling  maximum,  respec- 
tively. 

The  difference  exhibited  between  the  curves  may  be  due 
to  changes  in  the  character  of  the  air  mass  brought  about  by 
seasonal  migrations  of  the  quasi-permanent  thermal  low  over 
Nevada  and  the  North  Pacific  sub-tropical  high  located  west 
of  the  coast. 

Figure  5  was  used  to  calculate  a  revised  incident  radia- 
tion index  by  entering  with  the  newly-developed  upwelling 
index  values  and  by  applying  the  corresponding  percent  to 
the  theoretical  mean  monthly  radiation  for  clear  sky  condi- 
tions (possible  sunlight).   Figure  6  shows  the  possible  sun- 
light throughout  the  year  (from  Tont,  197  5)  as  Q  ,  and  the 

theoretical  radiation  at  the  surface  (Q.)  after  the  effects 

l 

of  fog  and  clouds  are  considered. 
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Figure  5.   Percent  of  possible  sunlight  as  a  function 
of  upwelling  index. 
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The  incident  solar  radiation  is  an  essential  part  of 
the  ecosystem  simulation  model  in  that  it  is  used  in  the 
equations  governing  the  production  of  organic  carbon  by  the 
photosynthetic  activity  of  phytoplankton.   These  equations 
are  discussed  in  detail  by  Pearson  (1975). 

C.   ALGAL  SINKING  TERM 

The  phytoplankton  sinking  term  which  was  added  to  the 
Pearson  (1975)  model  is  written  in  the  form  used  by  Riley 

(1965)  in  the  units  of  the  quantity  of  phytoplankton  trans- 

2 

ferred  per  day  (g  C/m  -day).   The  amount  of  phytoplankton 

that  sinks  out  of  the  mixed  layer  in  a  day  is  given  by: 

SINK  =  (V/Z)X2  (26) 

where:   SINK  =  flux  of  phytoplankton  out  of  mixed  layer 

(gC/m2-day) 

V  =  sinking  rate  (m/day) 

Z  =  mixed  layer  depth  (m) 

2 
X2  =  phytoplankton  biomass  in  mixed  layer  (gC/m  ) 

V/Z  =  fraction  of  phytoplankton  which  sinks  out  of 
mixed  layer  per  day 

Vertical  circulation  with  the  water  column  becomes 

significant  as  the  upwelling  current  speed  and  sinking  rates 

of  phytoplankton  approach  the  same  order  of  magnitude.  Since 

the  phytoplankton  are  non-mobile  and  are  generally  of  the 

same  density  as  the  water  or  have  developed  shapes  which 

increase  their  sinking  resistance,  they  will  be  carried  along 

with  vertical  currents  in  the  water  column.   The  vertical 

circulation  during  upwelling  opposes  sinking  but  may 
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accelerate  downward  transfer  during  brief  periods  of  surface 
convergence  and  associated  downwelling. 

The  complete  phytoplankton  sinking  term  is: 

SINK  =  ((V  -  W)/Z)X2  (27) 

where:       W  =  upwelling  speed  (m/day) 

An  average  value  of  one  meter  per  day  was  used "for  the  sink- 
ing rate  (V)  in  the  computer  simulation  model.   This  value 
was  determined  from  estimates  presented  by  Lehman  et  al . 
(1975)  and  Bannister  (1974).   The  SINK  term  acts  to  decrease 
the  phytoplankton  concentration  in  the  modeled  region  and  is 
subtracted  in  the  phytoplankton  equation  (2)  as  shown: 


^-rP-  =  PROD  -  RESP2  -  GRAZ  -  SINK  (28) 

dt 


D.   PHYTOPLANKTON  ADVECTION  TERM 

An  equation  describing  the  horizontal  advection  of  phyto- 
plankton is  included  in  the  revised  model.   This  term  is  in 
the  units  of  flux  and  is  written  as : 

ADVEC2  =  W  (X2)K  (29) 


where:     ADVEC2  =  the  change  of  concentration  of  phytoplank- 
ton over  time  (in  the  mixed  layer)  due  to 
advection  (g  C/m^-day) 

W  =  upwelling  speed  (m/day) 

X2  =  phytoplankton  concentration  in  the  mixed 
layer  (g  C/mO 

K  =  advection  coefficient  (m   ) 
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The  simplified  advection  term  used  in  the  model  is  de- 
rived from  the  equation  (21)  shown  in  the  background  theory- 
section,  which  is: 

Advective  change  =  W  (X2)  [-  ^  +  -^]  (21) 

The  Ekman  current  offshore  which  results  in  an  upwelling 
current  (W)  is  effective  over  the  depth  of  frictional  resist- 
ance (D) ,  but  the  phytoplankton  in  the  model  (which  are 
uniformly  distributed  in  the  mixed  layer)  will  be  advected 
at  the  average  velocity  over  the  depth  of  the  mixed  layer  as 
determined  by  the  fraction  of  the  mixed  layer  that  is  coinci- 
dent with  the  Ekman  layer.   It  is  known  that  the  mixed  layer 
varies  seasonally  from  about  10  to  100  meters  in  depth  and 
this  causes  an  annual  (seasonal)  variation  in  the  average 
horizontal  current  in  the  phytoplankton  advection  term.   The 
average  current  of  the  mixed  layer  can  be  determined  from: 

z 

U  =  —  /  U  dZ     where:  U?  =  a  function  of  e 

(30) 

A  shallow  mixed  layer  experiences  higher  average  veloci- 
ties than  a  deep  layer  as  seen  from  the  integral  expression 
(eq.  30). 

A  coefficient  to  describe  the  seasonal  average  of  the 
vertically  averaged  horizontal  current  in  the  mixed  layer  is 
included  in  the  advection  equation  as  shown: 

Advective  change  =  W  (X2)  [-  j^  Km  +  ^]         (31) 
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The  coefficient,  K  ,  represents  the  seasonal  average  of  the 
fraction  of  the  mixed  layer  which  is  coincident  with  the 
Ekman  layer  and  is  affected  by  the  Ekman  current. 

The  constants  in  the  advection  expression  (within  the 
brackets  of  eq.  31.)  are  lumped  into  the  advection  coefficient, 

K,  for  convenience  in  equation  (29).   The  range  of  K  varies 

-1       -2   -1 
from  about  10    to  10    m   depending  on  the  estimated  values 

of  the  gradient  distance,  (AX)  and  (AZ)  and  the  velocity  co- 
efficient, K  .   A  value  of  0.1  m   was  used  in  this  model 
'   m 

for  K. 


E.   PREDATION  LOSS 

The  predation  pressure  function  hypothesized  for  Monterey 
Bay  is  based  on  the  assumption  that  the  system  supports  a 
resident  population  of  herbivorous  zooplankton  predators 
throughout  the  year.   The  pressure  exerted  on  the  herbivore 
prey  by  this  maintenance  level  (N)  of  predators  follows  the 
curve  in  Figure  7  below  the  prey  density  threshold,  d.   The 
pressure  is  sufficiently  low  to  allow  growth  of  the  zooplank- 
ton stock,  but  high  enough  to  stabilize  growth.   A  second 
pressure  is  imposed  on  the  zooplankton  (see  Figure  7)  after 
their  density  reaches  the  critical  level,  d.   The  simulation 
approximates  conditions  which  may  be  brought  about  as  tran- 
sient predators  move  into  the  area  presumably  in  response  to 
increased  food  availability.   In  summary,  predation  pressure 
may  be  represented  by  a  two  part  function  with  a  pressure- 
density  curve,  a,  due  to  resident  predator  species  and  curve 
b   due  to  the  transient  predator  population. 
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Figure  7.   Predation  pressure  on  zooplankton  by 
transient  predators. 
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Predator  biomass  may  be  expressed  as  shown  in  Figure  8, 
as  a  simple  step  increase  or  some  other  function  of  prey 
density.   Consideration  of  the  step  type  function  is  impor- 
tant because  it  allows  the  predator  population  (and  conse- 
quent pressure)  to  maximize  during  those  seasons  when  the 
zooplankton  "standing  crop"  will  support  additional  numbers 
of  predators;  the  high  predator  pressure  will  deplete  the 
food  stock  to  a  lower  level  than  would  be  reached  by  the 
resident  predator  pressure  alone.   Measurements  of  zooplank- 
ton biomass  for  1974  (Traganza,  1975)  suggest  a  rapid  decline 
in  "standing  crop"  following  an  early  summer  maximum.   One 
might  suspect  that  a  transient  predator  population  is  forced 
out  of  the  region  once  the  food  sources  are  depleted.  This 
is  done  in  the  model  by  decreasing  the  predator  population 
when  the  zooplankton  biomass  declines  to  a  specified  level. 
The  predator  levels  were  related  directly  to  zooplankton 
"standing  crop"  by  triggering  an  increase  in  predator  popula- 
tion at  a  threshold  of  zooplankton  biomass  during  periods 
when  this  prey  population  was  on  the  upswing.   The  predators 
were  similarly  reduced  at  a  second  threshold  during  the  de- 
clining phase  of  zooplankton  growth.   A  set  of  four  condi- 
tional statements  is  used  in  the  simulation  to  describe  the 
predator-prey  relationship.   These  statements  specify  the 
following : 

IP   ^1   >   o    AND   X3  <  1.0   THEN:   FISH  =  1.0 
dt 

__   dX3_   >   0    AND   X3  >  1.0   THEN:   FISH  =5.0 

dt 
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Figure  3.   Predator  population. 
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IF      =$£      <       0         AND      X3    >    0.2      THEN:       FISH    =5.0 
at     —  — 


Ayr* 

IF   ^=~-   <   0    AND   X3  <  0.2   THEN:   FISH  =  1.0 
at  — 


The  modified  predation  pressure  term  in  the  zooplankton 
equation  of  The  model  is: 

LOSS  =  (L)(FISH)(X3)  (32) 


where 


LOSS  =  the  amount  of  herbivorous  zooplankton  biomass 
lost  per  day  as  a  result  of  predation 
(g  C/m2  day) 

L  =  loss  rate  or  percent  loss  per  FISH  per  day 
(day-1) 

FISH  =  number  of  predators  (non-dimensional) 

2 
X3  =  zooplankton  "standing  crop"  (g  C/m  ) 
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IV.   RESULTS 

A.   UPWELLING  CALCULATIONS 

The  mesoscale  wind  data  used  to  calculate  u'pwelling  for 
this  simulation  is  shown  by  Figure  9.   The  plot  depicts  the 
frequency  of  occurrence  (f  =  N/SN)  for  4-5-degree  segments 
of  the  compass.   The  prominent  wind  direction  is  as  usual 
from  the  northwest  .during  the  April  to  September  upwelling 
period  and  shifts  to  the  south  in  the  first  and  last  quarters , 
but  this  year  there  was  a  significantly  high  frequency  of 
northwesterly  winds  during  the  last  months  of  the  year.  The 
computer  generated  upwelling  index  (Figure  10)  shows  a  pri- 
mary upwelling  maximum  in  May,  a  second  peak  in  about  mid 
September  and  another  in  mid  November.   The  possibility  of 
a  secondary  upwelling  peak  is  suggested  by  temperature  obser- 
vations by  Traganza  et  al .  (1976).   Figure  11  shows  a  rise 
in  the  isotherms  peaking  30  to  45  days  after  the  indicated 
wind  initiated  upwelling  maximums  in  May  and  September.   Al- 
though the  water  column  does  not  respond  instantaneously  to 
surface  wind  stress,  the  delay  noted  here  is  excessive 
(Barham,  1957).   The  delay  may  be  attributable  to  the  assump- 
tion that  the  winds  measured  at  Monterey  Airport  are  actually 
representative  of  the  winds  over  the  bay,  when  in  fact  they 
may  not  be.   The  delay  may  also  be  caused  by  the  lack  of  suf- 
ficient data  to  accurately  depict  the  seasonal  trends  of 
temperature  in  the  water  column. 
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Traganza  et  al.  (1976)  observed  a  sharp  rise  in  phosphate 
(Figure  12)  with  the  major  upwelling  in  May  and  a  slight  in- 
crease in  phosphate  and  salinity  (Figure  13)  in  September. 
No  salinity  data  were  available  for  the  first  half  of  the 
year.   There  was  no  clear-cut  correlation  of  either  phosphate 
or  salinity  to  the  upwelling  index  during  November.   Nutrient 
data  from  four  to  nine  stations  taken  during  seven  cruises 
in  1974  were  averaged  (see  Figure  12)  and  show  a  rapid  rise 
in  October  which  is  probably  not  related  to  upwelling.  From 
the  five  year  study  of  Monterey  Bay  by  Bolin  (19  64)  it  is 
known  that  the  nutrient  concentration  is  characteristically 
low  over  the  depth  of  the  euphotic  zone  (0-200  meters)  for 
two  or  three  months  at  the  end  of  the  year.   The  restoring 
of  the  nutrient  level  in  all  but  the  upper  20  to  30  meters 
to  half  of  its  May  upwelling  value  while  the  surface  tempera- 
ture reached  apeak,  may  have  signalled  the  beginning  of  the 
Davidson  current  period  (see  Bolin  and  Abbott,  1962).   The 
Davidson  current  brings  a  southerly  winter  oceanic  water 
mass  into  the  Monterey  region.   This  "Davidson  water"  is 
characterized  by  lower  surface  temperatures  and  surface  salin- 
ities (due  to  high  amounts  of  rain) .   The  deeper  water  of 
the  euphotic  zone  has  higher  salinities  and  nutrient  concen- 
trations than  exist  at  Monterey.   The  mechanism  which  estab- 
lishes these  characteristics  may  be  winter  storm  mixing  oc- 
curring south  of  Monterey  or  upwelling  and  mixing  from  below 
brought  about  by  divergent  (cyclonic)  eddies  formed  in  the 
current  stream  as  it  moves  northward  along  the  coast 
(Traganza  et  al.,  1976). 
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Figure  13.   Seasonal  salinity  variation  off  Monterey 
Bay.   Values  are  the  mean  of  four  to  nine 
stations  (Traganza  et  al. ,  1976). 
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Traganza's  (19  76)  data  were  obtained  at  stations  in  an 
X  pattern  referenced  to  a  drogue  (Figure  14- )  and  at  inter- 
vals of  about  two  miles,  two  and  one-half  hours  apart.   The 
data  from  four  to  nine  stations  were  averaged  to  obtain  the 
values  shown  in  Figures  11,  12  and  13. 

Barham  (19  57)  argues  that  the  vertical  circulation  of 
Monterey  Bay  during  the  upwelling  period  is  characterized 
by  a  region  of  surface  divergence  coincident  with  the  head 
of  the  Monterey  submarine  canyon.   The  effect  is  due  to 
"channeling"  of  upwelled  replacement  water  along  the  canyon 
axis  with  a  plume  appearing  at  the  surface  near  shore  (see 
Figure  15).   Evidence  in  support  of  this  circulation  concept 
is  given  by  a  distinct  gradient  of  phytoplankton  concentra- 
tion in  the  surface  waters  outlying  the  canyon  head  (Barham, 
1957).   Barham  believes  that  the  gradient  is  due  to  the 
rapid  horizontal  advection  of  phytoplankton  innoculum  repre- 
senting a  potential  bloom  away  from  the  canyon  head  before 
the  population  shows  any  significant  growth.   During  periods 
of  reduced  upwelling  as  indicated  by  the  trends  of  the  iso- 
therms ,  phytoplankton  counts  over  the  canyon  rose  to  higher 
values ,  indicating  that  surface  divergence  or  horizontal  ad- 
vection had  diminished.   This  possibility  should  not  alter 
the  timing  of  the  model  which  is  related  to  the  general  up- 
welling, but  it  does  suggest  serious  spatial  sampling  prob- 
lems . 
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Figure  14.   Chart  of  the  Monterey  area  showing  station 
location. 
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Figure  15. 


Possible  "fountainhead"  effect  over 
canyon  head. 
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B.   SENSITIVITY  ANALYSIS 

The  model  was  tested  under  various  conditions  of  phyto- 
plankton sinking  rates,  predation  pressure  on  herbivorous 
zooplankton,  and  advection  of  phytoplankton  to  determine  the 
effect  on  the  simulation  of  the  seasonal  cycle  phosphate, 
phytoplankton  and  herbivorous  zooplankton. 

Losses  of  phytoplankton  by  sinking  and  therefore  food 
limitation  on  zooplankton  was  examined  by  setting  the  preda- 
tion terms  on  zooplankton  to  zero  and  varying  the  rate  of 
sinking  of  algal  cells,  i.e.  the  availability  of  food.   Sink- 
ing rates  of  zero,  three,  and  six  meters  per  day  were  used 
(Parsons  and  Takahashi,  1973  and  Riley,  1965).   The  effects 
are  shown  in  Figures  15,  17,  and  18.   Under  conditions  of 
zero  sinking,  a  rapid  rise  in  zooplankton  "standing  crop" 
in  late  summer  reduces  the  phytoplankton  biomass  quickly  to 
a  minimum  on  day  215.   The  large  decrease  in  phytoplankton 
allows  the  nutrient  concentration  to  remain  at  a  relatively 
high  level  (Figure  16)  since  the  uptake  of  nutrients  by 
phytoplankton  is  decreased  while  ongoing  zooplankton  excre- 
tion and  mixing  by  upwelling  add  to  the  nutrient  concentra- 
tion.  The  rate  of  growth  of  the  herbivorous  zooplankton 
population  or  "standing  crop"  is  shown  to  decrease  with  pro- 
gressively higher  phytoplankton  sinking  rates  (Figure  18) 
as  does  the  magnitude  of  the  maximum  zooplankton  biomass 
reached  during  the  year.   The  maximum  also  occurs  later  with 
increased  sinking  rate. 
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The  sinking  rate  was  varied  in  a  second  test  but  preda- 
tor pressure  was  exerted  on  the  zooplankton  by  employing 
the  step  function  predator-prey  relationship  detailed  earli- 
er.  Additional  predators  were  introduced  in  the  simulation 

o 
at  a  zooplankton  biomass  level  of  1.0  g  C/m   and  were 

removed  when  zooplankton  biomass  fell  to  0.2  g  C/m  .   Thresh- 
old values  were  determined  empirically.   The  zooplankton  bio- 

2  2 

mass  maximum  is  reduced  from  1-10  g  C/m   to  0.1-1.0  g  C/m 

by  the  addition  of  predation  pressure  (Figure  21).   The 
phosphate  curves  (Figure  19)  show  a  marked  response  to  phyto- 
plankton  growth  as  evidenced  by  a  lower  nutrient  minimum 
following  the  peak  of  the  summer  bloom.   A  considerable 
change  in  the  zooplankton  response,  i.e.  lower  and  earlier, 
is  evidently  due  to  reduction  of  food  sources  as  phytoplank- 
ton  is  allowed  to  sink  out  of  the  mixed  layer.   The  initial 
effect  of  increasing  the  sinking  .rate  from  zero  to  three 
m/day  on  zooplankton  is  a  shift  in  the  occurrence  of  the  peak 
to  approximately  70  days  later.   Further  increase  of  the 
sinking  rate  to  six  m/day  results  in  a  twofold  decrease  in 
the  carbon  biomass  of  zooplankton  with  an  additional  2  0  day 
delay  of  the  maximum.   It  is  apparent  that  the  rate  of  growth 
of  zooplankton  is  slowed  down  and  the  maximum  biomass  occurs 
later  because  of  food  limitation  but  that  predator  pressure 
is  responsible  for  limiting  the  magnitude  of  the  zooplankton 
biomass.   The  zooplankton  peaks  also  occur  earlier  when  pre- 
dation is  applied  and  also  occur  before  the  phytoplankton 
peaks  (Figure  20).   The  simulation,  therefore  requires  a 
positive  sinking  rate. 
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The  advection  rate  of  phytoplankton  was  next  varied 
while  the  sinking  rate  was  held  constant  at  1.0  m/day 
(Bannister,  1974)  and  predator  pressure  was  imposed  in  the 
same  manner  (and  with  the  same  thresholds)  as  in  the  second 
test.   The  results  are  shown  in  Figures  22,  23,  and  24.   A 
slower  rate  of  growth  of  zooplankton  is  produced  again  when 
the  advection  coefficient,  K,  is  increased  from  0.0  to  0.3 
m   .   Increasing  the  value  of  K  has  the  direct  effect  of  in- 
creasing the  rate  of  phytoplankton  advection.   The  single 
maximum  of  zooplankton  simulated  under  conditions  of  high 
advection  where  K  =  0.3  m    is  due  apparently  to  the  sensi- 
tivity of  the  model's  zooplankton  growth  equation  to  the 
mixed  layer  temperature  and  which  permits  a  rapid  growth 
of  zooplankton  during  the  mixed  layer  temperature  minimum. 
The  single  peak  of  zooplankton  (Figure  24)  coincided  with  a 
relative  temperature  minimum  about  day  2  60.   In  Pearson's 
simulation  a  ten  percent  decrease  in  mixed  layer  temperature 
had  a  marked  effect  on  zooplankton  growth.   This  effect 
carried  over  to  the  partially  modified  simulation  model. 
The  effect  of  rapid  zooplankton  growth  in  a  low  temperature 
mixed  layer,  despite  contradicting  /actors  of  growth  such 
as  low  food  availability,  is  due  to  the  sensitivity  of  the 
respiration  term  in  the  zooplankton  growth  equation  to  tem- 
perature.  This  condition  is  a  peculiarity  of  the  model  and 
it  is  doubtful  whether  the  real  ecosystem  behaves  in  the 
same  manner.   The  zooplankton  peak  precedes  the  phytoplankton 
peak  when  advection  is  zero  indicating  the  model  needs  an 
advection  term. 
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The  predation  pressure  terms  in  the  zooplankton  equation 
were  varied  (Figures  25,  26,  and  27)  while  keeping  the  sink- 
ing rate  at  1.0  m/day  and  the  advection  coefficient  at  K  = 
0.1m   .   The  conditions  of  the  test  were: 

1.  predation  pressure  set  to  zero; 

2.  predation  pressure  defined  by  the  linear  function  of 
zooplankton  biomass:  predator  pressure  =  L  times 
zooplankton  biomass,  where  L  =  0.01,  and  the  units 
are   g  C/m2-day  =  (day-1)(g  C/m2); 

3.  the  linear  function  in  2,  but  L  =  0.02; 

4-.   the  step  function  discussed  under  the  methods  section 
of  this  thesis. 

In  all  cases ,  the  resulting  zooplankton  growth  appears  food- 
limited  until  about  day  12  0  when  phytoplankton  growth  begins 
to  show  a  rate  change.   Complete  removal  of  predator  pressure 
allows  the  zooplankton  to  increase  rapidly  until  food  limita- 
tion (be  depletion  of  phytoplankton)  again  occurs  on  day  275. 

When  the  predator  pressure  increases  linearly  (with  L  =  0.01), 

2 
a  single  zooplankton  peak  of  3.05  g  C/m   occurs  about  day 

300.   This  is  not  consistent  with  observed  conditions 
(Traganza,  1976).   Doubling  the  predation  pressure  rate  to 
L  =  0.02  day"   severely  restricts  zooplankton  growth  but  pro- 
duces a  curve  with  a  hint  of  temporal  conformity  to  the 
observed  zooplankton  maximum  and  minimum  biomass  levels. 
When  the  predator  function  is  triggered  at  thresholds  of 
zooplankton  biomass,  the  zooplankton  exhibit  a  double  peak 
which  is  suggested  by  Traganza' s  (1976)  data. 
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C.   COMPARISON  OF  SIMULATION  RESULTS  WITH  OBSERVED  DATA 

The  computer  simulation  was  run  with  the  new  upwelling 
and  radiation  indices  and  the  step  function  predation  pres- 
sure term  detailed  in  the  preceding  section.   The  sinking 
rate  was  set  at  1.6  m/day  and  the  advection  coefficient  set 
to  K  =  0.1  m"1. 

A  comparison  of  the  simulation  results  with  data  observed 
by  Traganza  et  al .  (1976)  is  shown  in  Figures  26,  27,  and  28. 
The  calculated  nutrient  .values  lag  the  averaged  mixed  layer 
concentrations  by  20  to  30  days  during  the  summer  months  with 
a  moderate  amount  of  error  in  magnitude.   The  nutrient  maxi- 
mum occurs  on  day  144  with  a  value  of  2.08  yg-at  P/l  in  the 
simulation.   A  rapid  decrease  in  phosphate  levels  coincident 
with  the  summer  phytoplankton  bloom  displays  the  impact  of 
nutrient  uptake  by  phytoplankton  in  the  model. 

The  simulated  phytoplankton  peak  occurs  40  days  after 
the  nutrient  peak  during  the  period  of  maximum  incident  radia- 
tion (see  Appendix  C).   The  modeled  response  of  phytoplankton 
was  relatively  inflexible  with  regard  to  the  radiation  index, 
i.e.  the  phytoplankton  peak  was  found  to  occur  within  a  few 
days  of  the  time  of  maximum  incident  radiation  when  the  sink- 
ing rate  was  varied  from  about  1.0  m/day  to  1.8  m/day.   The 
relationship  of  the  phytoplankton  peak  to  the  radiation  peak 
is  explained  by  examination  of  the  phytoplankton  growth  equa- 
tions (shown  in  Appendix  B).  A  slight  increase  in  phytoplank- 
ton biomass  was  simulated  during  December  (day  330)  which 
could  be  due  to  the  increased  phosphate  concentration  at  the 

end  of  the  year. 
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Twenty  days  after  the  phytoplankton  maximum,  the  model 

simulated  a  zooplankton  peak  with  a  concentration  of  0.97 

2 
g  C/m  .   The  zooplankton  biomass  growth  corresponds  well 

with  the  observed  data  until  approximately  the  beginning  of 

September.   The  measured  values  increase  to  a  maximum  of 

2 
1.85  g  C/m   in  mid-December  while  the  model  simulated  a 

2 
recovery  peak  of  O.M-7  g  C/m   and  slightly  earlier. 
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V.   DISCUSSION 

The  Monterey  upwelling  ecosystem  simulation  developed  in 
this'  thesis  has  been  shown  to  follow  the  seasonal  trends  of 
the  observed  phosphate,  phytoplankton,  and  zooplankton  data. 
There  are  errors  in  magnitude  of  the  simulated  response  at 
various  times  of  the  year.   The  differences  between  observed 
and  simulated  zooplankton  are  most  likely  due  to  the  fact 
that  the  model  simulates  herbivorous  zooplankton  while  the 
observations  include  both  herbivores  and  carnivores .   Quanti- 
tative data  defining  the  seasonal  ratio  of  herbivorous  to 
carnivorous  zooplankton  are  needed  to  verify  the  simulation 
results.   The  encouraging  aspect  of  the  model  is  its  ability 
to  follow  the  seasonal  trends  and  the  fact  that  the  timing 
of  the  response  of  phytoplankton  biomass  to  the  nutrient  con- 
centration and  the  response  of  zooplankton  to  the  phytoplank- 
ton cycle  is  biologically  sound,  i.e.  a  moderate  delay  is 
simulated  between  the  peaks  of  successive  trophic  levels. 

Since  all  the  factors  likely  to  affect  the  three  state 
variables  (phosphate,  phytoplankton  and  zooplankton)  are  not 
included  in  the  model,  some  error  in  magnitude  should  appear 
In  recalling  the  philosophy  expressed  earlier,  the  objective 
of  this  thesis  has  been  to  create  a  time  simulation  of  the 
dynamics  of  phosphate,  phytoplankton  and  herbivorous  zoo- 
plankton in  a  limited  area  over  a  long  time.   This  objective 
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has  been  met  by  balancing  the  four  model  characteristics  of 
generality,  realism,  precision  and  simplicity. 

From  studies  by  Barham  (19  57)  and  Bolin  (19  64)  it  is 
noted  that  there  may  be  characteristic  seasonal  patterns  of 
plankton  dynamics  in  response  to  characteristic  patterns 
of  the  hydrography  of  the  Monterey  upwelling  ecosystem.   It 
is  perhaps  important  to  realize  that  the  apparent  limita- 
tions of  the  model  as  developed  this  far  are  relatively 
small  when  the  use  of  the  simulation  to  reproduce  the 
characteristic  patterns  is  considered  per  s_e. 

As  more  and  better  data  with  which  to  define  the  dynamics 
are  obtained,  the  simulation  model  of  the  Monterey  upwelling 
ecosystem  can  progress  to  a  more  refined  level. 
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VI.   SUGGESTIONS  FOR  FURTHER  RESEARCH 

The  following  comments  are  made  as  a  result  of  questions 
which  arose  during  the  modification  of  the  Monterey  ecosystem 
model.   This  list  is  provided  to  identify  some  areas  which 
may  warrant  further  investigation.   It  is  by  no  means  com- 
plete, but  should  serve  as  a  guide. 

1.  The  most  significant  outcome  of  this  research  has 
been  the  recognition  of  the  importance  of  accurate  forcing 
functions.   The  presence  of  the  submarine  canyon  in  Monterey 
Bay  and  the  likely  effect  on  upwelling  suggests  further 
effort  be  directed  to  verifying  the  upwelling  index. 

2.  The  Monterey  Bay  area  is  not  homogeneous  in  physical, 
chemical  and  biological  parameters  as  evidenced  by  the 
patchiness  of ' biological  samples  and  the  observable  spatial 
gradients  of  the  physical  and  chemical  properties. 

3 .  Additional  work  on  the  biological  aspects  of  the 
model  including  verification  and  refinement  of  the  predation 
terms  and  development  of  an  appropriate  kinetic  expression 
is  indicated. 

4.  The  simulation  appears  overly-sensitive  to  tempera- 
ture variations  in  the  mixed  layer. 

5.  Various  investigators  have  shown  the  importance  of 
phytophagous   fish  in  ecosystem  models.   A  fish  herbivore 
term  should  be  investigated. 
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6.  Some  question  as  to  the  stability  of  the  CSMP  360 
routine  over  long  time  simulation  has  been  raised.   This 
might  be  investigated  and  resolved  by  initializing  the 
model  part  way  through  the  year. 

7.  The  effect  of  varying  the  forcing  functions  in  the 
revised  model  has  not  yet  been  studied. 

8.  The  advection  expression  might  be  expanded  to  in- 
clude the  effects  of  all  current  systems  influencing  the 
Monterey  region. 
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